function cfdAssembleIntoGlobalMatrixFaceFluxes
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018 
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   This function assembles algebraic equation coefficients from the
%   contribution of the face fluxes of the current term of the equation
%--------------------------------------------------------------------------

% Get fluxes
theFluxes = cfdGetFluxes;

% Get info
theNumberOfInteriorFaces = cfdGetNumberOfInteriorFaces;
theNumberOfFaces = cfdGetNumberOfFaces;
owners = cfdGetOwnersSubArrayForFaces;
neighbours = cfdGetNeighboursSubArrayForFaces;

upperAnbCoeffIndex = cfdGetUpperAnbCoeffIndex;
lowerAnbCoeffIndex = cfdGetLowerAnbCoeffIndex;

% Get coefficients
theCoefficients = cfdGetCoefficients;

ac = theCoefficients.ac;
anb = theCoefficients.anb;
bc = theCoefficients.bc;

% Assemble fluxes of interior faces
% 对于anb（a_neighbor）的解释：相当于以稀疏的形式来存储Ax=b的稀疏矩阵A，刨去了对角线
% 这是一个cell结构，有numberOfElements行，
% 每行是一个数组，以cell索引iEle为例，按照elementFaces(iEle)内的邻接面索引顺序
% 依次存储着与各个邻居element交互的系数（根据这两个element共享的面的通量）

for iFace=1:theNumberOfInteriorFaces % 遍历每个内部面
    own = owners(iFace);
    nei = neighbours(iFace);
    % 下面两行是执行从全局索引（iFace）到局部索引（分别在owner和neighbor中各一套局部索引）的转换
    iOwnerNeighbourCoef = upperAnbCoeffIndex(iFace);% iOwnerNeighbourCoef表明iFace这个（全局）索引的面在它的owner cell的面集合里是第几个面（局部的索引）
    iNeighbourOwnerCoef = lowerAnbCoeffIndex(iFace);% iOwnerNeighbourCoef表明iFace这个（全局）索引的面在它的neighbor cell的面集合里是第几个面（局部的索引）
    
    % Assemble fluxes for owner celliNeighbourOwnerCoef
    % 注意：对于owner，ac(own)是加上FluxCf，anb{own}(iOwnweNeighborCoeff)是加上FluxFf
    ac(own)                       = ac(own)                       + theFluxes.FluxCf(iFace);% 应该是对角线元素（系数）
    anb{own}(iOwnerNeighbourCoef) = anb{own}(iOwnerNeighbourCoef) + theFluxes.FluxFf(iFace);% 应该是表示除对角线外的稀疏矩阵
    bc(own)                       = bc(own)                       - theFluxes.FluxTf(iFace);% 应该是右端项？但为什么是FluxTf而不是FluxVf呢：因为uFVM求解的是残差形式的方程Ax=b-Ax*
    
    % Assemble fluxes for neighbour cell
    % 注意：对于neighbor，ac(nei)是减去FluxFf，anb{nei}(iNeighborOwnerCoeff)是减去FluxCf，与owner相反！
    ac(nei)                       = ac(nei)                       - theFluxes.FluxFf(iFace);
    anb{nei}(iNeighbourOwnerCoef) = anb{nei}(iNeighbourOwnerCoef) - theFluxes.FluxCf(iFace);
    bc(nei)                       = bc(nei)                       + theFluxes.FluxTf(iFace);
end

% 关于方程右端项bc为什么是用FluxTf而不是FluxVf，是因为与书上不同，uFVM求解的方程不是Ax=b
% 而是对应的残差方程：A(x^{n} + x^{*}) = b，其中x^{n}是上一迭代步的结果，待求的是增量x^{*}
% 也即是当前解x^{n}和最终稳态收敛解之间的残差，因此形式上为：Ax^{*} = b - Ax^{n}，
% 此处的bc也就是b - Ax^{n}。注意到b包含原PDE右端源项（案例中可能是0）和离散后的非正交项（-FluxVf），
% 而Ax^{n}也就用FluxCf(iFaces).*T(owners_f)+FluxFf(iFaces).*T(neighbours_f)表示，
% 对应参照cfdAssemble_xxx_DiffusionTerm，所以这里是bc -= FluxTf
% FluxTf 表示穿过这个面的整个通量

% Assemble fluxes of cfdBoundary faces 边界面只有一个owner element需要管
for iBFace=theNumberOfInteriorFaces+1:theNumberOfFaces
    own = owners(iBFace);
    
    % Assemble fluxes for owner cell
    ac(own) = ac(own) + theFluxes.FluxCf(iBFace);
    bc(own) = bc(own) - theFluxes.FluxTf(iBFace);
end

% Store
theCoefficients.ac = ac;
theCoefficients.anb = anb;
theCoefficients.bc = bc;

cfdSetCoefficients(theCoefficients);

end
